Collapse of triaxial bright solitons in atomic Bose-Einstein condensates 
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We study triaxial bright solitons made of attractive Bose-condensed atoms characterized by the 
absence of confinement in the longitudinal axial direction but trapped by an anisotropic harmonic 
potential in the transverse plane. By numerically solving the three-dimensional Gross-Pitaevskii 
equation we investigate the effect of the transverse trap anisotropy on the critical interaction strength 
above which there is the collapse of the condensate. The comparison with previous predictions [Phys. 
Rev. A 66, 043619 (2002)] shows significant differences for large anisotropies. 

PACS numbers: 03.75.Lm,03.75.Kk,03.75.Hh 



The experimental achievement of quantum degeneracy with ultracold alkali-metal atoms [l|, Q has opened the 
possibility of studying various topological configurations of the Bose-Einstein condensate (BEC) with repulsive or 
attractive inter-atomic interaction Q . Dark solitons in repulsive BECs have been experimentally achieved ten years 
ago [H, while bright solitons have been detected only more recently in two different experiments @] involving 
attractive BECs of 7 Li vapors. In these latter experiments, an optical red-detuned laser beam generated along the 
axial direction of the sample is used to trap the attractive BEC by a cylindric isotropic transverse confinement; the 
BEC propagates along the longitudinal axis of the cylinder without relevant spreadings. Recently, also 85 Rb atoms 
have been used to achieve the Bose-Einstein condensation and investigate the formation of bright matter-wave solitons 
during the collapse 0- 

Many theoretical works have been devoted to the study o f cig ar-s hap ed and axially symmetric bright-soliton config- 
urations, also in presence of an axial periodic potential EE EI El El EI El EE E3 • The transverse confinement 
produced by the isotropic harmonic potential in the cylindric radial direction plays a crucial role in giving rise to single 
1, I s , 10, TlL 12 1 or multiple (l3l, [T3| metastable bright solitons, which collapse above a critical number of particles 
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12L Il4|. These theoretical investigations showed that increasing the inter-atomic strength, e.g. by Fano- 
Feshbach resonances, makes the bright soliton less cigar-shaped. In particular, a quasi-spherical shape is achieved 
only when the interaction strength approaches the critical value that signs the collapse [lCJ, E2l ■ 

In this paper we study an attractive BEC trapped by an anisotropic harmonic potential in the transverse plane 
and without confinement in the axial direction. Under this trapping condition the BEC admits stable bright-soliton 
configurations, which are are generally triaxial. The deformation of the transverse anisotropic confinement can be 
described by two independent parameters [l8j or by a unique quantity, the ellipticity (l9j . In both cases, by using 
the numerical integration of the three-dimensional Gross-Pitaevskii equation, we investigate as a function of the 
transverse trap anisotropy the critical interaction strength above which the triaxial soliton collapses, i.e. it shrinks 
to the zero-size ground-state of infinite negative energy. We compare our results with previous numerical predictions 
EH and find that our stability domain is significantly smaller. 

Let us consider an attractive BEC without confinement in the axial direction z and confined in the transverse plane 
[x, y) by the anisotropic harmonic potential 

U(x,y) = ^(^+^y 2 ) , (1) 

where m is the mass of a Bose-condensed atom, and uji, ll>2 are the two frequencies of the harmonic confinement. 
With the aim of working in scaled units we set 

LJl = Ai Wj_ , lo 2 — A 2 wj_ . (2) 

In particular, if aj_ = {%/ (muj^)) 1 / 2 is used as characteristic harmonic length of the system, then lengths may be 
measured in units of a± and energies in units of Hlu±. 

The dynamics of an attractive BEC can be accurately described by the adimensional time-dependent 3D Gross- 
Pitaevskii equation (3D GPE), given by 
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where ^(r, t) is the macroscopic wave function of the condensate and 
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FIG. 1: Integrated density profiles p(x), p(y) and p(z) of the triaxial bright soliton along the three Cartesian axes obtained by 
the numerical integration of the 3D GPE. Interaction strength g — 2N\a s \/a± = 1.2 and trap ellipticity e = 0.4. 



is the interaction strength, with TV the number of atoms and a s < the s-wavc scattering length of the inter-atomic 
potential. Setting 



*(r,t) = V(r) e - *"* , 



from Eq. ((3|) one finds the stationary 3D GPE 



(5) 



ip = fi ip , (6) 



(7) 

Stable solutions of Eq. © are called bright solitons [1, [H, [l(], EL E3 ■ They correspond to an attractive BEC with a 
self-confinement along the z axis. 

To determine the solutions of Eq. (J5]) we solve Eq. by using a finite-difference Crank-Nicolson algorithm with 
imaginary time [2l[ and a spatial mesh of 160 x 160 x 160 points (for details see the Appendix). In this way we 
determine the wave function ^>(r) of the metastablc bright soliton and we can also calculate the integrated density 
profiles of the triaxial bright soliton along the three spatial directions, given by 



where the chemical potential pi is fixed by the normalization 

|V^(r)| 2 d 3 r = 1 



p(x) - 

p{y) : 
P{z) 



|7/>(r)| 2 dy dz , 
IVK 1 ")! 2 dx dz , 
^(r)! 2 dx dy 



(8) 
(9) 
(10) 



To study the critical strength above which there is the c ollap se, we consider first the simpler case of elliptic 
transverse confinement. As explained by Jamaludin et al. [19{ . it is possible to consider an elliptic transverse 
harmonic confinement and consequently to parametrize the transverse anisotropy of the harmonic confining potential 
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FIG. 2: Critical strength g c for the collapse of the triaxial bright soliton as a function of the ellipticity e of the elliptic transverse 
harmonic confinement, Eq. (jXTJ . Filled circles: numerical results obtained with the 3D GPE. Dashed line: prediction of Eq. 

G3. 
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FIG. 3: Chemical potential fi as a function of the interaction strength g obtained with the 3D GPE. e = 0.8 is the ellipticity 
of the elliptic transverse harmonic confinement. 



(fTj) by using a unique parameter, the trap ellipticity e. In terms of the ellipticity e, the scaled harmonic frequencies 
are written as 

\! = VT~e, X 2 =VT+e, (11) 

with e restricted to the interval [—1,1]. Clearly e = corresponds to the isotropic transverse confinement, while 
e = ±1 implies the absence of confinement along the x axis (e = 1) or along the y axis (e = —1). 

As an example, in Fig. [T]we plot the density profiles of the bright soliton choosing the interaction-strength g = 1.2 
and an elliptic transverse confinement with ellipticity e = 0.4. The shape of the bright soliton strongly depends on 
the ellipticity e of the transverse potential and the interaction strength g. By varying e and g the bright soliton can 
be spherical-shaped, cigar-shaped, disk-shaped, but also fully triaxial. 

Our numerical investigation shows that under the condition e > the width a x of the soliton along the x axis is 
always close to 1 (in units of a±). The width a y of the soliton along the y axis is equal to a x only for e = 0; moreover 
<j y becomes extremely large as e — > 1. Obviously, with e < the behaviors of a x and a y are interchanged. The width 
<j z along the z axis is instead controlled by the interaction strength g: a small g implies a very large <j z , while when 
g is sufficiently large the width a z is around 1. In addition it exists a critical strength g c above which there is no 
solution, i.e. the wave function of the metastable soliton collapses to the zero-size ground-state of infinite negative 
energy. 

In Fig. [2] we plot this critical strength g c as a function of the ellipticity e of the transverse trap. Our numerical 
results based on the integration of the 3D GPE arc displayed as filled circles. The figure shows that when the trap 
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FIG. 4: Critical strength g c as a function of the scaled frequency A2 with Ai = 1. The interaction strength is g — 2iV|a s |/oj 
Filled circles: numerical results obtained with the 3D GPE. Dashed line: prediction of Eq. (|12p . 



is perfectly symmetric (e = 0) the critical strength g c reaches its minimum value, g c = 1.35. Instead, as |e| — > 1~ 
the critical strength has its maximum value given by g c = 1.66. We notice that when |e| — > 1 _ , the frequency of 
confinement along one of the two transverse directions goes to zero, but only at e = ±1 the triaxial bright soliton 
becomes unbounded. 

It is interesting to compare our results with previous predictions based on numerical calculations and scaling [Til ] . 
According to these predictions [Til ] the critical strength g c is simply given by the formula 
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In Fig. [2] the dashed line is obtained with Eq. (fl~2"]) and the scaled frequencies given by Eq. (fTTj) . The figure shows 
that there are relevant differences between our numerical results (filled circles) and Eq. (fl~2|) for |e| close to 1. In fact, 
Eq. (fl"2"]) implies that g c — > +00 for |e| — > l - . Actually, Eq. (fl"2"|) is based on the hypothesis of a attractive BEC with 
triaxial harmonic confinement of frequencies Ai, A2 and A3 under the conditions Ai,A2 3> A3 and A3 — > [TT|] , but 
these conditions on the harmonic frequencies break down for |e| close to 1. As previously, discussed, our numerical 
results suggest instead a finite value of g c for |e| — > 1 _ , as confirmed also by a fully Gaussian variational approach 
0. 

An important issue is the dynamical stability of the triaxial bright solitons we have found. According to the 
Vachitov-Kolokolov criterion j22j, the fundamental solitons are stable if they satisfy the condition dfi/dg < 0. We 
have verified that up to the collapse this condition is always satisfied by our bright solitons. For completeness, in Fig. 
[3] we show the calculated chemical potential (i versus the interaction strength g for the triaxial bright solitons with 
ellipticity e = 0.8. 

Let us now investigate the general case where Ai and A2 are independent. Keeping fixed one of the two harmonic 
frequencies, e.g. Ax, we may independently tune the other, A2. Without loss of generality we fix Ai = 1. We find that 
the critical strength g c approaches a maximum finite value when the trapping frequency A2 tends to zero. This effect 
is shown in Fig. [4J where we plot the critical strength g c as a function of A2 with Ai = 1. Instead, for large values 
of A2 the critical strength g c becomes smaller. By using a Gaussian variational approach [l8| we have indeed verified 
that g c — ► as A2 — > +00. For the sake of completeness, in Fig. 4 we have included also the prediction of Eq. JT2]) 
with Ai = 1. Remarkably, there are deviations not only for small values of A2 but also for large values of A2. 

In conclusion, in this work we have investigated the collapse of triaxial bright solitons in Bose-condensed atoms 
under transverse anisotropic harmonic confinement by using the 3D Gross-Pitaevskii equation. Our predictions on the 
stability domain of these triaxial bright solitons can be useful for future experimental investigations with deformed 
atomic waveguides. 

This work has been partially supported by Fondazionc CARIPARO. The authors thank Boris Malomcd and Flavio 
Toigo for useful suggestions. 



5 



Appendix: Numerical method 



The numerical integration of the time dependent GPE, in Eq. ((3]), is obtained by using a finite-difference Crank- 
Nicolson scheme modified with a split operator technique, adapted to the integration of a Schrodinger equation [2lj . 
This approach has been successfully applied in various problems and configurations [Til , [l5| . 

First, we write Eq. ([3]) in the form 

ift£*(r,t) = (H^t) + H 2 (r,t) + H 3 (r,t))V(r,t), (13) 



dt 



where 



H a (r,t) ee + U(x a ) - i. 9 |*(r,t)| 2 , (14) 

with a = 1,2, 3 and xi = x, X2 = y, 2:3 = 2. Here we have used the fact that the external potential is separable: 
U(x,y,z) = U(x) + U(y) + U(z). In this way we split the full Hamiltonian in three sub-Hamiltonians, so that at 
each time we have to write the Laplacian with respect to one coordinate only, leading to the solution of a tridiagonal 
system, and to huge savings in computer memory. 

Equation (fTB")) is integrated using the splitted Crank-Nicolson scheme 

*(r,t + ft) = _L_(l-A l( t)/2)x 

(l-A 3 (t)/2)x 



l + A 3 (t)/2 
1 

l + Ai(t)/2 



(l-A a (t)/2)*(r,t). (15) 



where S t is the integration time step, and A a (t) = iStH a (r, t)/Ti. The splitting is carried out so that the commutators 
are exact up to the order S\ included. There is obviously a problem with the nonlinear term ^^(r, t)| 2 , because we 
should really use a '5 somehow averaged over the time step St, not a '3/ evaluated at the beginning of the time step. 
To circumvent this problem, we use a predictor-corrector step, where each integration step is really done in two times: 
going from the time t to the time t + S tl the first time we used ^(r, t) in the nonlinear term, obtaining a "predicted" 

*P(r, t + St); we then repeated the integration step, starting again from ^(r, t), but using | ^(r, t) + ^(r, t + 5)j in 

the nonlinear term. 

In our numerical method the wave function is discretized in the following way 

*(r,t) = ${x i ,y j ,z k ,t s ) (16) 

where x % = xq + i 5 X , y = yo + j S y , z k = zq + k S z , and t s — s S t , with i, j, k, s integer numbers. Second derivatives 
are approximated by finite-difference formulas. For instance, along the x axis we use 

5 2 , j k ^^ *(x*+ 1 ,y3,z k ,t s )-2*(x\y>,z k ,t s ) + *(x l ~ 1 ,yi,z k ,t s ) 
dx 2 ( ' r ' ' >~ 81 ' [ ' 

In the upper panel of Fig. [5] we show the typical behavior of the energy E of the triaxial bright soliton as a function 
of the imaginary time t. We use a triaxial Gaussian as initial trial wave function ^(r, t = 0) = ^initiaii^), normalizing 
to one the norm of ^(r, t) at each time step. The energy E of the system, given by 



E = J **(r,t) 



iv 2 + i(A 2 , 2 + A^ 2 )-i(2 7 r g )|*(r,t)| 2 



*(r,i)d 3 r, (18) 



decreases during the (imaginary-)time evolution and eventually reaches its asymptotic value Efi na i. We find that the 
asymptotic value Efi na i depends of the number N s x N s x N s of points in the spatial mesh. Nevertheless, as shown 
in the lower panel of Fig. for sufficiently large values of iV s the energy Efi na i saturates to the exact value. As a 
final check, we have verified that smaller values of the time step St, with respect to the one we use (St = 0.05), do not 
modify the final results within the third digit. Moreover, we have checked that the final energy E^ na x and the final 
wave function ^ final (j) do not depend on the initial trial wavefunction ^ initial (r) . ^finaiij) is the wave function 
of triaxial bright soliton. In the case of collapse, we find that the final energy is Efi na i = — 00 and the final wave 
function is a Dirac delta peak, centered at r = 0. 
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FIG. 5: Triaxial bright soliton with g = 1.2 and e = 0.8. Upper panel: energy E of the bright soliton as a function of the 
imaginary time t for 3 values of N s . The spatial mesh has N 3 x N s x N a points. Lower panel: asymptotic energy Ef ina i of the 
bright soliton as a function of N s . 

A very recent and complete review of numerical methods used to solve the Gross-Pitaevskii equation has been 
written by Muruganandam and Adhikari [23| . In this paper the finite-difference Crank-Nicolson scheme we have used 
is explained with many details. 
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